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We study the accretion on to the neutron star in Be /X-ray binaries , using a 3D 
SPH code and the data imported from a simulation by lOkazaki et al.l l)2002|) for a 
coplanar system with a short period (Porb = 24.3 d) and a moderate eccentricity 
(e = 0.34), which targeted the Be/X-ray binary 4U 0115+63. For simphcity, we adopt 
the polytropic equation of state. We find that a time-dependent accretion disc is formed 
around the neutron star regardless of the simulation parameters. In the long term, the 
disc evolves via a two-stage process, which consists of the initial developing stage and 
the later developed stage. The developed disc is nearly Keplerian. In the short term, 
the disc structure modulates with the orbital phase. The disc shrinks at the periastron 
passage of the Be star and restores its radius afterwards. The accretion rate on to the 
neutron star is also phase dependent, but its peak is broader and much lower than 
that of the mass-transfer rate from the Be disc, unless the polytropic exponent is as 
large as 5/3. Our simulations show that the truncated Be disk model for Be/X-ray 
binaries is consistent with the observed X-ray behaviour of 4U 0115+63. 

Key words: accretion, accretion discs - hydrodynamics - methods: numerical - 
binaries: general - stars: emission-line. Be - X-rays: binaries 



1 INTRODUCTION 

The Be/X-ray binaries represent the largest subclass of high- 
mass X-ray binaries. About two-thirds of the identified sys- 
tems fall into this category. These systems consist of, gen- 
erally, a neutron star and a Be star with a cool (~ 10^ if) 
equatorial disc, which is geometrically thin and nearly Ke- 
plerian. The orbit is wi de (16 d < Porb 5, 243 d) and usually 
eccentric (e>0.3) (e.g., |zioiw23|2TO^. 

Most of the Be/X-ray binaries show only transient ac- 
tivity in the X-ray emission and are termed Be/X-ray tran- 
sients. Be/X-ray transients show periodical (Type I) out- 
bursts, which are separated by the orbital period and have 
the lumiocity of Lx ~ 10'^®~'^^ergs~^, and giant (Type 
II) outbursts of Lx ^ lO'^^erg s~^ with no orbital modula- 
tion. The transient nature in the X-ray activity of Be/X- 
ray binaries is considered to result from the interaction 
between the accreted material from the circumsteller disc 
of the Be star and the rotat ing magnetized neutron star 
JStella. White fc RosneJllQSeH . If the accreted material is 
dense enough to make the magnetospheric radius smaller 
than the corotation radius, the accretion on to the neutron 
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star causes a bright X-ray emission (the direct accretion 
regime). Otherwise, the magnetospheric radius is larger than 
the corotation radius, and the accretion is prevented by the 
centrifugal inhibition. The system is then in quiescence (the 
propeller regime). 

During the outbursts, the neutron stars in most Be/X- 
ray binaries show spin-up (e.g., 4U 115+63. GS 18 43-02 
EXO 2030+375 and XTE J1946-f-26) jBildsten et al. 199 
iFineer et all 1 19991: IWilson et ai]l2002l: IWilson et al. .200 
and some even show quasi-periodic o scillations (QPOs 
(e.g. , A 0535+ 26 and EXO 2030+375) iFineer et alJ Il9& 
Wils on et al.l 12002") . These observational features strongly 
suggest the presence of an accretion disc during the out- 
bursts. 

The accretion disc is also indicated to exis t between the 
outbu rsts in seve ral systems (e.g.. A 0535+26) jFinger et alJ 
[199^ see a lso [ Ziolkowski l2002tl . On the other hand, 
Iciark et alJ Jl99ar searched the optical/infrared contribu- 
tion from the accretion disc in A 0535+26 during the X-ray 
quiescent phase, and found no signature of the accretion disc 
in propeller regime. Obviously, more work is needed to deter- 
mine whether the accretion disc also exists in the quiescent 
state. 

Recently, lOkazaki et alJ (l2002l) studied the interaction 
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Figure 1. Illustrative diagram of a Be/X-ray binary. At the pe- 
riastron passage, the gas in an outer part of the Be-star disc is 
transferred to the neutron star and forms an accretion di sc. Th e 
right panel is taken from a simulation by Okaza ki et all J2002f) . 
while the left panel is from one of our simulations. 



between the Be-star disc and the neutron star, using a three- 
dimensional (3D) smoothed particle hydrodyna mics (SPH) 
code iBenzlll990l : iBate. Bonnell fc Pried Il995ll . They also 
qualitatively discussed the relation between the X-ray lu- 
minosity and the accretion rate using the simulated mass- 
capture rate by the neutron star. In their simulations, how- 
ever, the neutro n star was modelled by a sink particle 
jBate et al.lll995^ with the size of the Roche lobe. Hence, 
no information was available how the captured material ac- 
cretes on to the neutron star. Moreover, no direct compari- 
son with the observed X-ray data was possible. 

In this paper, we study the accretion on to the neu- 
tron star in Be/X-ray binaries, using a 3D SPH code. The 
simulation should ideally take account of all the processes 
at work in a Be/X-ray binary, such as the disc evolution 
around the Be star, the mass transfer from the Be-star disc 
to the neutron star, and the accretion on to the neutron 
star. Such a simulation, however, would require very intense 
computations. 

Therefore, we confine ourselves to simulate only the ac- 
cretion flow around the neutron star, with an outer bound- 
ary condition constructed by the data of the particles cap- 
tured by the neutron star from a high-resolution simula- 



tion by lOkazaki et al.l l|2003) • This situation is illustrated 
in Fig. 1, where the right panel taken from their simula- 
tion shows the mass transfer from the Be-star disc to the 
neutron star at periastron, while the left panel shows the 
accretion disc formed around the neutron star in one of our 
simulations. 



2 NUMERICAL MODEL 
2.1 SPH code 

Simulations presented here were performed with a 3D SPH 
co de. The SPH c ode i s basically the same as that used 
bv lOkazaki et alJ i2002f) . w hich is based on a version orig - 
inally developed by Benz jBend Il99d : iBenz et alJ 11990) . 
The smoothing length is variable in time and space. The 
SPH equation with the standard cubic-spline kernel are inte- 
grated using a second-order Runge-Kutta-Fehlb erg integra- 
tor w ith individual time steps for each particle iBate et alJ 
1^9^, which results in an enormous computational saving 
when a large range of dynamical time-scales are involved. 



In our code, the accretion flow on to the neutron star is 
modelled by an ensemble of gas particles, each of which has 
an arbitrary and negligible mass chosen to be 10~ ^^Mq with 
a variable smoothing length. Note that the accretion flow is 
nonself-gravitating in our simulations. On the other hand, 
the Be star a nd the neutr on star are modelled by two sink 
particles ( Bat e et al]ll995l) with corresponding masses. The 
gas particles which fall within a specifled accretion radius are 
accreted by the sink particles. We assume that the neutron 
star has the fixed accretion radius of 3.0 — 5.0 x 10~^a, a 
being the semi-major axis of the binary. For the Be star, 
we adopt a variable accretion radius of 0.8rL, where tl is 
the Roche-lobe radius for a circular binary. This is because 
our interest is only in the fiow around the neutron star. An 
approximate formula for tl is given by 



rL 



: 0.462 



1 + q 



1/3 



(1) 



(see, for example. IWarnerlll995l : lFrank. King fc Raineir2002h 

with the mass ratio q = Mk/M*, where Mx and M* are the 
masses of the neutron star and the Be star, respectively, and 
D is the distance between the stars. 

In disc simulations using a three dimensional SPH code 
with the standard cubic spline kernel, one has an approx- 
imate relation between the Shakura-Sunyaev viscosity pa- 
rameter ass and the linear SPH artificial viscosity parame- 
ter QfsPH as 



ass 



1 h 



(2) 



if the non-linear SPH artifical viscosity parameter /3sph is 
set equal to zero (Oka zaki ot a, l. 20Q2j. Here h is the SPH 
smoothing length and H is the scale-height of the disc. 

In a simulation we report in this paper, we adopt 
aspH ~ 1 and /3sph = 2 in which case ass is variable in 
time and space. In the other simulations, however, we adopt 
a constant value of ass in order to roughly model the a 
disc. In these simulations, qsph = lOassH /h was variable 
in time and space and /3sph ~ 0. 

In our simulations, we adopt a polytropic equation of 
state specified by the exponent F, in which the effects of an 
energy dissipation resulting from viscosity and a radiative 
cooling are approximately boxed. 



2.2 Boundary condition 

Our binary configuration is the same as that of 
j^kazaki ct al. (2002). We set the binary orbit on the x-y 
plane with the semi-major axis along the a;-axis. Initially, the 
neutron star is at the periastron. It orbits about the Be star 
with the orbital period Porb = 24.3 d and the eccentricity 
e = 0.34. The unit of time adopted here and throughout the 
paper is Porb. The Be disc is coplanar with the binary orbital 
plane. As the Be star, we take a BOV star with M, = ISMq, 
P* = 8Rq and T^s = 26000 K, where R, and Tefi are the 
stellar radius and effective temprature, respectively. For the 
neutron star, we take Mx = 1.4Mq and Px = 10® cm. 
These parameters are adopted to target 4U 0115-1-63, which 
has shown only Type II X-ray ou tbursts followed by a sh ort 
series of Type I X-ray outbursts (iNegueruela et al]l2001^ ■ 

As mentioned earlier, we concentrate only on the sim- 
ulations of the accretion flow on to the neutron star. This 
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Orbital Phase 

Figure 2. Orbital-phase dependence of tlie mass-capture rate 
by tlie neutron star^ taken from a high-resolution simulation by 
lOkazaki et al.l j20o3)- The data is folded on the orbital period 
oyer 32 < t < 47. The mass-capture rate is measured in units of 
P_iiMq yr~^, where p— ii is highest local density in the Be-star 
disc normalized by 10~^^gcm~^, a typical value for Be stars. The 
right axis shows the X-ray luminosity corresponding to the mass- 
capture rate with the X-ray emission efficiency t] = I, where r] is 
defined by Lx = '?MxMcap/J?x- 

is done by constructing an outer boundary condition from 
the data of particles capt ured by the neutron st ar in a high- 
resolution simulation by lOkazaki et alJ ll2002h . which was 
run for < t < 47 and finally had ~ 140000 particles with 
ctss = 0.1 . We first fold the data on the orbital period 
over 32 < t < 47 to reduce the fluctuation noise. The phase 
dependence of the mass-capture rate by the neutron star 
obtained by this procedure is shown in Fig. |21 Note that 
there is a strong phase-dependence in the mass-transfer rate. 
The peak comes slightly after the periastron passage and 
decreases by two orders of magnitude by apastron passage. 
We then fit the spatial and velocity distribution of the cap- 
tured particles in each phase bin with separate Gaussian 
functions. This is for running simulations with an arbitrary 
mass-transfer rate from the Be disc. We finally model the 
mass-transfer process from the Be disc to the neutron star 
by injecting gas particles at a given phase-dependent rate 
shown in Fig. 2 with the spatial and velocity distributions 
specified by these Gaussian functions. The injected particles 
are assumed to have the initial temperature of Tcff/2. 

As described earlier, the mass of each gas particle is 
arbitrary and negligible in our simulations. All the mass- 
related quantities discussed later, such as the disc mass and 
the mass- accretion rate, are normalised by using the mass- 
capture rate shown in Fig. 2. 



2.3 Possibility of an accretion disc formation 

Before describing the results of our simulations, we discuss 
below the possibility of the formation of an accretion disc 
around the neutron star. 

We first consider the radius -Rcirc at which the mate- 
rial transferred from the Be disc is circularized via the in- 
teraction with the previously transferred gas. We assume 
the circularization radius -Rcirc to be approximately given 
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Orbital Phase 

Figure 3. Orbital-phase dependence of the circularization radius 
(the upper panel) and the ratio of the viscous time-scale to the 
orbital period (the lower panel) for the gas transferred from the Be 
dis c in 4U 0115-1-63, base d on the data imported from a simulation 
by Okaz aki et al.l j20o2) . The horizontal dotted line in the upper 
panel denotes the corotation radius Rq = 6.0 X 10~^a. 

by -Rcirc = J^/GMx, where J is the initial specific angular 
momentum of the transferred material. 

In the upper panel of Fig.|3 we show the dependence of 
^icirc on the orbital phase, based o n the data taken from a 
simulation bv lOkazaki et all (|20o3). It is noted that there is 
a strong phase-dependence of 7?circ . It has a minimum of ~ 
lO^^a slightly before periastron and then rapidly increases 
with phase. For comparison, we also show the corotation 
radius, i?n, for 4U 0115+63 by the dotted line. Here, Rq is 
given by 

where Pspin = 3.61 s is the spin period of the neutron star in 
4U 0115+63. Note that Rdrc > Rn for most of the orbital 
phase. Therefore, we expect that an accretion disc is formed 
around the neutron star in 4U 0115+63. 

Next, we discuss whether the accretion disc formed 
around the neutron star is persistent or transient. This is 
done by comparing the viscous accretion time-scale with the 
orbital period. For simplicity, we assume the accretion disc 
to be geometrically thin and isothermal with the Shakura- 
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Table 1. Summary of models of our simulations. The initial number of particles is 301 and the mean mass injection rate is 2.5 X 
10~^^p— llM0yr~^ in all simulations. The first column represents the model. The second column is the polytropic exponent T and the 
third column is the number of SPH particles at the end of the run, which is given in units of Parh in the fourth column. The viscosity 
parameters adopted are given in the fifth column and the sixth column is the inner radius of the simulation region. The seventh column 
is the number of injected particles for one orbital period. The last two columns are the number of accreted particles for the last one 
orbital period and the corresponding mean accretion rate, respectively. 



Model 


Polytropic exponent 


A^SPH 


Run time 


Viscosity parameters 


Inner boundary 










r 


(final) 


(Porb) 




'"in/a 






(p_liM0yr-i) 


1 


1.2 


45687 


8 


ass = 0.1 


3.0 X 10"^ 


8657 


2990 


8.3 X 10-^^ 


2 


1.2 


35821 


8 


ass = 0.1 


5.0 X 10-3 


8717 


4054 


1.1 X 10-" 


3 


1 (isothermal) 


39295 


8 


ass = 0.1 


5.0 X 10-3 


8844 


2676 


7.4 X 10-12 


4 


5/3 (adiabatic) 


20118 


8 


ass = 0.1 


5.0 X 10-3 


8594 


6851 


1.9 X 10-" 


5 


1.2 


47730 


8 


asPH = l,/3sPH = 2 


5.0 X 10-3 


8705 


2671 


7.5 X 10-12 


6 


1.2 


17488 


1 


ass = 0.1 


5.0 X 10-3 


34778 


16829 


1.2 X 10-" 



Sunyaev viscosity parameter ass and evaluate the accretion 
time-scale at r = i?circ. The ratio of these two time-scales 
for 4U 0115+63 is then given by 



-forb 



1 



27rQssCs 



1.2 X 10' 



i / QSS \ ~ 

Vo.i / 



1+i 

1/2 



1/2 



10*K 



(4) 



where Cs is the sound speed and Td is the disc temprature. 
The orbital-phase dependence of r^is/ PoA is shown in the 
lower panel of Fig. |H] It is immediately observed that the 
viscous accretion time-scale is much longer than the orbital 
period for most of the orbital phase. Therefore, once formed, 
the accretion disc around the neutron star in 4U 0115-1-63 
will survive over the whole orbital phase. Note that this 
result depends little on the orbital period, because i?circ in- 
creases almost linearly with increasing a. 



3 ACCRETION FLOW AROUND THE 
NEUTRON STAR 

Based on the procedure described in the previous section, 
we have performed several simulations with different pa- 
rameters. Details of the models are given in Table Q The 
polytropic exponents considered are F = 1 (model 3), 1.2 
(models 1, 2, 5 and 6) and 5/3 (model 4). We have taken the 
inner radius of the simulation region, rin, to be 5.0 x 10~^a 
in all simulations but model 1 with a smaller inner radius of 
Hn = 3.0 X IQ-^a, which was run to study the effect of the 
inner radius on the simulation results. In models 1-4 and 6, 
we have assumed ass = 0.1, where aspH varies in time and 
space and /3sph = 0. For comparison purpose, we have also 
run model 5 with qsph = 1 and /3sph = 2, where ass varies 
in time and space. Model 6 was run solely for clearly illus- 
trating the initial build-up phase of an accretion disc. For 
this purpose, we have taken the injection rate of the SPH 
particles in model 6 four times as high as in other models 
and run only for one orbital period. 

Analysing multi-wavel ength long term monitor ing ob- 
servations of 4U 0115-1-63, iNeeueruela et al.l J200J) found 
that the Be star undergoes quasi-cyclic (~ 3 — 5yr) activity, 
losing and reforming its circumsteller disc. Then, the mass 
transfer to the neutron star should vary, depending on the 
state of the Be disc. By the time the Be disc is fully de- 
veloped and becomes capable of supplying mass to the neu- 
tron star, even a large accretion disc would be exhausted. 



Therefore, to investigate the evolution of an accretion flow 
around the neutron star after the Be disc is fully developed, 
we performed 3D SPH simulations based on the numerical 
models described above, where there is intially no gas par- 
ticle around the neutron star. Since all models have shown 
similar structure and evolution of the accretion flow, we will 
mainly show the results from model 1 with F = 1.2 and 
nn = 3 X lO'^a. 



3.1 Initial phase of the disc formation 

We consider the initial phase of accretion flow around the 
neutron star. As discussed in Section r2.3l it is expected that 
the material transferred from the Be disc near the periastron 
does not directly fall on to the neutron star and forms a 
disc-like structure with a size corresponding to its specific 
angular momentum. 

For illustration purpose, we give in Fig. |l|the snapshots 
of the initial accretion flow around the neutron star from 
model 6, which covers t = Q — 1. At t = Q, the Be star is 
at periastron. Each panel shows the surface density in the 
range of four orders of magnitude in the logarithmic scale. 
We note that the accreting matter first forms an elliptical 
ring around the neutron star, which is then relaxed to an 
elliptical disc via the viscous diffusion. Because the viscous 
time-scale is longer than the orbital period, the orbits of the 
accreting particles are not circularized and the flow is highly 
eccentric during the initial phase of accretion. Comparing 
the number of particles, A'^sph, with the number of particles 
inside the effective Roche radius of the neutron star, A^Rocho, 
we also note that most of the material transferred from the 
Be disc stays inside the effective Roche radius of the neutron 
star. It should be noted that these dynamical features are 
common to all models studied in this paper. 

Fig-Elshows the time-dependence of the radial structure 
of the accretion flow around the neutron star at the same 
times as in Fig. 2] In each panel, the solid line, the dashed 
line and the dash-dotted line denote the surface density E, 
the azimuthal velocity normalized by the Keplerian ve- 
locity at the inner boundary and the radial Mach number 
Vr/cs, respectively. In the figure, the density is integrated 
vertically and averaged azimuthally, while the velocity com- 
ponents are averaged vertically and azimuthally. From Fig.|5] 
we note that while the azimuthal velocity distribution grad- 
ually approaches the Keplerian one, the behaviour of the 
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Figure 4. A sequence of snapshots showing the initial phase = — 1) of the disc formation around the neutron star in model 6 with 
r = 1.2 and 

'"in — 5 X 10 '^CL. The origin is set on the neutron star. Each panel shows the surface density in a range of four orders 
of magnitude in the logarithmic scale. The dashed circle denotes the effective Roche lobe of the neutron star. The short line-segment 
indicates the direction of the Be star. The unit of time is the orbital period Porb- Annotated at the top-left corner of each panel are the 
number of SPH particles, AfspHi and the number of particles inside the effective Roche lobe of the neutron star, A'^Rochc- 




r/a r/a r/a r/a 



Figure 5. Time-dependence of the radial structure of the accretion flow shown in Fig.|4] In each panel, the solid line, the dashed line 
and the dash-dotted line denote the surface density S in units of p_iigcm~^, the azimuthal velocity normalized by the Keplerian 
velocity at the inner boundary and the radial Mach number Vr/c^, respectively. For comparison, the Keplerian velocity distribution is 
shown by the dash-three-dotted line. The density is integrated vertically and averaged azimuthally, while the velocity components are 
averaged vertically and azimuthally. 



radial velocity is very complex, reflecting the dynamical pro- 
cess of the disc formation shown in Fig. 2] 

The surface density distribution also shows the initial 
evolution of the accretion flow described above. In the panel 
at t = 0.12, there appears a hump at r ~ 0.08a caused by 
the infall of injected particles. Outside the hump the density 
decreases rapidly, making a ring-like structure with a sharp 



outer edge. This rapid decrease in the density is due to the 
ram pressure of the infalling particles. As the hump moves 
inward, the density gradient outside the hump becomes less 
steep. By the next periastron passage {t = 1), a disc-like 
structure is formed, the inner radius of which is at r ~ 0.02a. 

In our simulations, the presence of the inner simulation 
boundary causes an artificial decrease in the density near the 



6 Kimitake Hayasaki and Atsuo T. Okazaki 



boundary, because the gas particles which pass through the 
boundary are removed from the simulation. As a result, the 
density distribution has a break near the inner boundary, 
as seen in Fig. |3 We have found that our simulations are 
reliable for > 2rin. 



3.2 Structure of a developed disc 

Many Be/X-ray binaries with known orbital parameters 
have orbital eccentricities in the range of 0.3 to 0.9. In 
such systems, the mass-transfer rate from th e Be disc is ex- 
pecte d to have the strong phase dependence ijOkazaki et alJ 
l2002h . Accordingly, the structure of an accretion disc formed 
around the neutron star in these systems is also likely to be 
strongly phase dependent. In this section, we first discuss the 
phase-dependent structure of a developed accretion disc in 
model 1 and then compare that with those in other models. 



3.2.1 Model 1 

Fig. |S| gives a sequence of snapshots of the accretion flow 
around the neutron star covering one orbital period for 7 < 
t < 8. The format of the figure is the same as that of Fig. |4| 
except that the surface density is shown in the range of three 
orders of magnitude. We note from the figure that the disc 
shrinks at periastron passage by a negative torque exerted by 
the Be star and restores its radius afterwards by the viscous 
diffusion. We also note that the disc is significantly eccentric 
even after it is fully developed, although the eccentricity 
is much smaller than that in the initial phase of the disc 
formation. 

Fig. |7| shows the radial disc structure at the orbital 
phases corresponding to Fig. |^ The format of the figure 
is the same as that of Fig. 5. It is noted from Fig.|7|that the 
disc is nearly Keplerian at any phase and the accretion flow 
is highly subsonic except in the outermost part for a short 
period of time when the material transferred from the Be 
disc falls towards the neutron star. We call this state of the 
accretion disc developed. We have found that these features 
of a developed disc hold for t > 5 in model 1. As in the initial 
phase of the disc formation, the ram pressure by the super- 
sonic infall of matter transferred from the Be disc around 
periastron enhances the density in the outermost part of 
the accretion disc, making the disc outer edge sharp (see 
panels for 7.12 <t< 7.37 in Figs. |H] and 0. Afterwards, the 
disc gradually expands until the next periastron passage of 
the Be star. 

Fig. |H| shows how hot the disc is at t = 7.5 in models 1 
and 3-5. In each panel, the solid line, the dash-dotted line 
and the dashed line denote the disc temprature normalized 
by the initial temprature To — 1.3 x 10* K, the ratio of 
the smoothing length to the harf-thickness of the disc h/H 
and the relative disc thickness H/r, respectively. From the 
upper- left panel of Fig.|Hl it is noted that the disc in model 1 
is vertically resolved except for the innermost region of the 
disc (r<5.0 x 10~^a). The disc temprature varies from ^ 
1Q*K at ~ O.lo to ~ 6.5 x 10'' K at the inner boundary. 
The disc is geometrically thin with the relative thickness of 
about 0.1-0.2. 



3.2.2 Models 2-5 

In this section, we describe the effects of the simulation pa- 
rameters in Table on the structure of the accretion flow 
around the neutron star. 

We ran model 2 to study the effect of the inner bound- 
ary on the disc strcture. Fig.|n|shows the snapshots and the 
radial structures of the accretion fiow around the neutron 
star at t = 7 (periastron) and t = 7.5 (apastron) for model 2, 
which has the same simulation parameters as model 1, ex- 
cept for a larger inner boundary of rin = 5.0 x 10~^. Com- 
paring Fig. |U] with the panels at the corresponding times in 
Figs. |S| and |7| for model 1, we note that the presence of the 
inner boundary artificially reduces the density in a narrow 
region (r < 2rin) near the boundary, but outside of it, the 
disc structure is similar for both models. 

In order to study the dependence of the disc structure 
on the polytropic exponent, we also ran models 3 (F = 1; 
isothermal case) and 4 (F = 5/3; adiabatic case), for which 
the polytropic exponent of 1.2 for models 1, 2 and 6 is in 
between. Fig. 1101 shows the snapshots and the radial struc- 
tures of the accretion disc for model 3 (isothermal case) at 
the same times as in Fig. |U] From the figure, we note that 
the distributions of the surface density and both of the az- 
imuthal and radial velocity components are similar to those 
in the initial phase of the disc formation for F = 1.2 (see 
Figs. 2]and|21l. Consequently, in model 3 the accretion disc 
is still developing at t ~ 7. This is because the viscous time- 
scale of the isothermal disc, particularly in an inner region, 
is longer than those in other models (see also Fig. I^J. From 
the upper-right panel of Fig. |S| it is noted that the disc in 
model 3 is not vertically resolved in most of the disc region. 

Fig. II II shows the disc structure for model 4 (adiabatic 
case). We note that the disc is nearly Keplerian at both 
orbital phases and that the surface density is lower and the 
disc is bigger than those for model 2 with F = 1.2. This is due 
to a higher viscous stress in model 4, resulting from a higher 
pressure. As expected, in model 4, the accretion rate is much 
higher than the models with lower polytropic exponents, 
which we will see in a later section. From the lower-left panel 
of Fig. |5] it is noted that the disc in model 4 is much hotter 
than in other models, reaching the temprature of ~ 5.0 x 
10® K near the inner boundary. The disc is vertically resolved 
except for the innermost region. 

We also studied a model with the standard, artificial 
viscosity parameters. Model 5 has the same simulation pa- 
rameters as those in model 2, except that it has constant 
articficial viscosity parameters aspH = 1 and /Ssph = 2, 
for which ass varies in time and space. Fig. 1121 shows the 
disc structure for model 5. From the figure and the lower- 
right panel of Fig. |H| it is immediately noted that the inner 
boundary affets the disc structure much more widely in this 
model. This is mainly because the inner disc region is not 
resolved vertically for this number of SPH particles. Since 
the ratio of the smoothing length to the harf-thickness of the 
disc h/H in the inner region is much higher than unity, the 
shear viscosity parameter ass in this model is much higher 
in the inner region than that in other models. 

Except for the above-mentioned dependence of the disc 
structure on the simulation parameters, the accretion discs 
in our models share the common features: They are nearly 
Keplerian after developed and their structures modulate 



Accretion disc formation in Be/X-ray binaries 7 




0.2 



0.1 r' 



N 'i^6154 t=7.Z5 



t 



-0.2 -0.1 0.1 0.2 
x/a 



0.1 



0.1 0.2 
/a 



-0.2 




-0.2 -0.1 0.1 0.2 
x/a 




Figure 6. A sequence of snapshots of a developed accretion disc in model 1, which cover one orbital period (t = 7 — 8). The format of 
the figure is the same as that of Fig.|l] except that the surface density is shown in a range of three orders of magnitude in the logarithmic 
scale. 
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Figure 7. Time-dependence of the radial structure of the developed accretion disc shown in Fig. [5] The format of the figure is the same 
as that of Fig. El 



with the orbital phase. They are also significantly eccen- 
tric. We will analyse the non-axisymmentry of the accretion 
discs in Be/X-ray binaries in a subsequent paper. 



is truncated via the resonant interaction with the neutron 
star, as long as ass <S 1, wh ich has recently been confirmed 
by numerical simulations bv lOkazaki et all (|20o3). 



5. .8. 5 Size of the accretion disc in Be/X-ray binaries 

lArtvmowicz fc Lubowl lll994f) investigated the 
tidal/resonant truncation of circumsteller and circumbinary 
discs in eccentric binaries and found that a gap is always 
opened between the disc and the binary orbit. Fo llow- 
ing their formulation^ Ne gueruela fc Okazakil ||200 j) and 
lOkazaki fc Neeueruel i 1200 J) have shown that the Be disc 



It is interesting to study whether the tidal/resonant 
torque by the Be star also truncates the accretion disc 
around the neutron star. In order to have a measure of 
the accretion disc radius at which the disc density has a 
major break, we have applied a non-linear least squares fit- 
ting method to the radial distribution of the azimuthally- 
averaged surface density S, adopting the following simple 
fitting function, 
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(r/rd)-^ 



1 + (r/rd)i 



(5) 



where p and q are constants and rd is the radius of ac- 
cretion disc. This m ethod is the same as that adopted by 
iQkazaki etaP l|20o3). 

Fig. ll3l shows the orbital- phase dependence of the accre- 
tion disc radius in model 1. To reduce the fluctuation noise, 
we folded the data on the orbital period over 5 < i < 8. It is 
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Figure 12. Same as Fig.|2l but for model 5. 



noted that the disc radius varies between 0.04a and 0.05a. 
The disc shirinks at periastron due to a negative torque by 
the Be star and gradually restores its radius towards apas- 
tron due to the viscous diffusion. The rapid increase of the 
disc radius before periastron is due to the infall of matter 
from the Be star. 

Next, we examine whether the disc radius shown in 
Fig. 13 is determined by the resonant truncation as is 
the case in the Be d i sc in Be/X-ray binaries. Fig. 8 of 
lArtvmowicz fc Lubowl (|l994^ gives the truncation radius of 
a circum-secondary disc as a function of the orbital eccen- 
tricity e and the disc Reynolds number Re for systems with 
the mass ratio of 0.1, which is similar to that of Be/X-ray 
binaries. Since the accretion disc in our model is geometri- 
cally thin and nearly Keplerian, the Reynolds number Re is 
given by Re = {l/ass)(r/ H)^ ~ (l/ass)(t'Kep/cs)^, where 
VKep is the Keplerian velocity. In model 1, Re ^ iq^-s.s 
for O.Ola < r < 0.1a. From Fig. 8 of Art vmowicz &: Lubowl 



10'' 



we have the trun- 



J1994) . with e = 0.34 and Re 
cation radius of ~ 0.14a, which corresponds to the 6:1 reso- 
nance radius. Note that this radius is much larger than the 
disc radius shown in Fig. 1131 It is likely that the interac- 
tion with the Be star is significantly non-linear owing to its 
large mass. Then, the accretion disc would be truncated well 
inside the 6:1 radius given above. Even such a strong inter- 
action, however, is unlikely to truncate the disc at a radius 
as small as that shown in Fig. 1131 

It is suggestive that the disc radius is about the same 
as the circularization radius -Rcirc of the matter transferred 
from the Be disc after the periastron passage (see Fig. I^J. 
Since the viscous time-scale at the circularization radius is 
much longer than the run time of our simulations, the mat- 
ter once circularized dose not have time enough to change 
its radius significantly. Therefore, we conclude that the disc 
radius shown in Fig. ll3l is determined by the specific angular 
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Figure 13. Orbital-phase dependence of the accretion disc radius 
r'tj. To reduce the fluctuation noise, the data was folded on the 
orbital period for 5 < t < 8. 




Figure 14. Evolution of (a) the disc mass and (b) the disc 
size normalized by the semi-major axis for < t < 8. In 
panel (a), the disc mass is measured in units of P-hIO^^^Mq. 



momentum of the matter transferred from the Be disc. It is 
probable, however, that the disc expands well beyond the 
circularization radius and is truncated by a tidal torque in 
a time-scale longer than the viscous one. 



3.3 Disc evolution 

As mentioned in Section 1, most of the Be/X-ray binaries 
show transient X-ray activities, of which nature is considered 
to result from the interactions between the ac creted mate- 
rial a nd the rotating magnetized neutron star iStella et alJ 
Il986l) . In this section, we first describe the growth of the ac- 
cretion disc in mass and radius. We then examine the evolu- 
tion of the mass- accretion rate obtained by our simulations. 
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from which we finally discuss the X-ray activity expected in 
Be/X-ray binaries. 

3.3.1 Growth of the disc 

Fig. 1141 al shows the evolution of the disk mass in model 1 
over the period of < i < 8. In the long term, the disc mass 
increases with time, because the the gas transferred from the 
Be star accumulates in the accretion disc owing to the vis- 
cous time-scale longer than the orbital period. On the other 
hand, in the short term, it modulates with the orbital phase. 
The disk mass starts increasing, when the mass-supply from 
the Be disc begins slightly before every periastron passage. 
It shows a rapid increase for a short while, followed by a 
gradual decrease by accretion with no mass-supply from the 
Be disc, which lasts until the next periastron passage. At 
first glance, the disc mass in model 1 seems to increase sec- 
ularly and approach no steady average state. On closer in- 
spection, however, we note that the increasing rate of the 
disc mass gradually decrases with time after the disc is de- 
veloped. This suggests an interesting possibility that the disc 
reachs a steady average state in a time-scale much longer 
than the term shown in Fig. 14(a). 

In Fig. I14r b'l. we show the evolution of the disc radius 
over < t < 8. One may expect that the disc radius gradu- 
ally increases with the disc mass. From the figure, however, 
we note that this is the case only for t>5. In fact, for t^5, 
the disc radius decreases as the disc grows. This counter- 
intuitive decrease in radius is related to the change in the 
disc eccentricity. In the early stage of disc formation, the 
disc/ring is highly eccentric (see Fig. 2J, for which the ra- 
dius rd calculated by equation Q is larger than that for 
a circular disc with the same disc mass. Since the disc ec- 
centricity decreases with time, so does the disk radius. In 
model 1, the effect of the disc circularization is balanced 
with that of the disc growth at t ~ 5. After this epoch, the 
disc radius increases with the disc mass. Thus, the evolution 
of the accretion fiow around the neutron star in Be/X-ray 
binaries involves a two-stage process, which consists of the 
initial developing stage and the later developed stage. 



The orbital modulation in the mass-accretion rate is 
significantly difi^erent between the developing disc and the 
developed disc. The first peak becomes weak as the disc 
developes. After the disc is fully developed, the second peak 
is higher than the first peak. 

Fig. llSf al also shows that, once developed, the accre- 
tion disc in 4U 0115+63 shows the peak X-ray luminosity 
of Lx ^ 2 ■ 10^^ ergs~^. This level of X-ray em ission enters 
the transition regime (see lCampana et al.ll20oil) and is con- 
sistent with the observed X-ray luminosity in the quiescent 
state. 

For comparison purpose, we present in Fig. llfcil the evolu- 
tion of the mass-accretion rate and the corresponding X-ray 
luminosity in models 2-4 for < t < 8. The thick and thin 
solid lines and the thin dotted line are for models 2, 3 and 4, 
respectively. We see from the figure that the accretion rate 
increases with increasing polytropic exponent. In models 2 
(r = 1.2) and 3 (F = 1), the accretion rate profile has the 
double peaks as described above. The long-term change in 
the peak profiles is also similar to that in model 1. We note 
that the disc in model 2 entered the stage of the developed 
disc at t ~ 6. On the other hand, the disc in model 3 is at 
the developing stage even at t = 8, which is consistent with 
the disc structure shown in Fig. IIUI 

In contrast to the models with low F, the accretin rate 
in model 4 has a single and periodic peak per orbit, even 
after the disc is developed. With F — 5/3, the pressure in 
the disc in model 4 is much higher than that in other models. 
This means that much larger viscous torque exerts on the 
disc in model 4 than in other models, because the viscous 
stress in the a-viscosity prescription is proportional to the 
pressure. As a result, the viscous torque dominates the tidal 
torque in model 4 even at periastron. This is why there is no 
peak at periastron and the accretion rate shows the regular 
orbital modulation with the amplitude much higher than in 
other models. 

It is important to note that the disc evolution and its 
orbital modulation is sensitive to the polytropic exponent 
F. This strongly implies that the cooling effciency in the 
accretion disc plays an important role in the X-ray activity 
in Be/X-ray binaries. 



3.3.2 Phase-dependent accretion 

Fig. I15f a'l shows the evolution of the mass-accretion rate 
and the corresponding X-ray luminosity in model 1 for < 
t < 8. Here we calculated the X-ray luminosity by Lx ~ 
rjGMxMacc/ Rx with the X-ray emission effciency rj — 1 for 
the neutron star. 

The thick and thin lines denote the mass-accretion rate 
in models 1 and 2, respectively. Note that the mass-accretion 
rate in our simulations generally has double peaks per orbit. 
The first peak, which occurs at the periastron, is due to the 
contraction of the disc by the tidal torque of the Be star, 
except that the peaks at periastron during the initial phase 
of disc formation are mainly due to the direct accretion of 
gas with low specific angular momentum. The second peak, 
which lags behind the peak in the mass-transfer rate from 
the Be disc shown in Fig. |5| results from the viscous accre- 
tion of matter transferred from the Be disc. The first peak 
could be artificial since it is related to the presence of the 
inner simulation boundary. 



3.3.3 Magnetospheric radius 

Changes in the mass flux on the magnetosphere of the neu- 
tron star determine whether the system is in the direct ac- 
cretion regime or in the propeller regime. While the matter 
falls directly on to the neutron star in the direct accretion 
regime, no material accretes on to the neutron star in the 
propeller regime by the centrifugal inhibitation. If the mag- 
netospheric radius _Rm is smaller than the corotation radius 
Rq, the system is in the direct accretion regime, otherwise, 
in the propeller regime. 

Assuming the neutron star magnetic field as the dipole- 
like field, the magnetospheric radius Rm, obtained by equat- 
ing the magnetic field pressure of the neutron star to the ram 
pressure of the accreting matter from the accretion disc, is 
written by 

\10i6 ergs-i J \MqJ 
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JPrinde fc Reesll972l:lLamb. Pethick fc PinesllOTl see also 
iFrank et alJl2002r i. where ^ is the magnetic moment. Here, 
A: is a constant model parameter of order unity, depend- 
ing on the physics and geometry of accretion. For spher- 
ical accretionj_fc_^ 0.91. The original model calculations 
by iGhosh fc LambI ^1979^ for disc accretion on to the 
aligned rotating magneti zed neutron star o btain ed k ~ 0.47, 
whereas calculations by lOstriker fc Shd lll995ll. which in- 
cluded an induced wind not considered bv l^^osh fc Larnbl 
il979tl . resulted in the value of fc ~ 0.923 for disc accretion. 

Fig. llSf b) shows the evolution of the magnetospheric 
radius -Rm normalized by the corotation radius Rn, where 
k = 0.47 was adopted. In the figure, the thick and thin 
lines are for models 1 and 2, respectively. Recall that these 
models have the same parameters except for the radius of the 
inner simulation boundary rin; Hn = 3 x 10^'^ in model 1, 
while rin = 5 X 10""^ in model 2. It should be noted that 
the magnetospheric radii in both models evolve in a very 
similar way, despite the different inner radii. Our result on 
the magnetospheric radius is robust in this sense. 

We also note from Fig. llSf b) that the ratio of the mag- 
netospheric radius to the corotation radius, Rm/Rsi, varies 
between ~ 2fc and ^4fc. Consequently, for a plausible range 
of fc (0.5 ^ fc ^ 1), the magnetospheric radius in models with 
r < 1.2 is almost always larger than the corotation radius 
so that the accretion on to the neutron star is inhibited, 
which is consistent with the obcrvod X-ray luminosity of 
4U 0115-1-63 in quiescence (Campana et al.. .2001j') . On the 
other hand, we can see from equation (O and Fig. 1161 that 
model 4 with F = 5/3 is consistent with the observed X- 
ray behaviour of 4U 0115-1-63 only if fc>0.6. For fc ~ 0.5, 
the neutron star in model 4 enters the direct accretion 
regime after every peristron passage of the Be star. Then, 
the system would show regular, periodic X-ray outbursts at 
a level of ~ 5 x lO'^^ergs"^, which has not been observed for 
4U0115-^-63. 



4 SUMMARY AND DISCUSSION 

For the first time, we have performed numerical simulations 
of accretion on to the neutron star in Be/X-ray binaries, tak- 
ing into account the phase-dependent mass transfer from the 
Be disc. We have adopte d the mass-transfer ra te from a high- 
resolution simulation bv lOkazaki et al] ll2002h for a coplanar 
system with a short period (Porb = 24.3 d) and moderate ec- 
centricity (e — 0.34), which targeted 4U 0115+63, one of the 
best studied Be/X-ray binaries. We have found that a nearly 
Keplerian, non-steady accretion disc is formed around the 
neutron star regardless of the parameters adopted. The for- 
mation of the accretion disc involves a two-stage process, 
which consists of the initial developing stage and the later 
developed stage. The disc structure shows a strong depen- 
dence on the orbital phase, because both the tidal potential 
and the mass-transfer rate change periodically. The disc is 
also significantly eccentric, because the gas particles from 
the Be disc originally have elliptical orbits. 

In 4U 0115-1-63, Type I X-ray outbursts have occurred 
only as a short series after oc casional Type II X-ray out- 
bursts (INegueruela et al.ll200ll) . Except for such occasions. 



the system is in the q uiescent state (10^^ — 10"^^ ergs~^) 
iCampana et alJ 1200 J) . The peak X-ray luminosity esti- 
mated from our simulations corresponds to the transition 
regime between the direct accretion regime and the propeller 
regime and is consistent with the observed X-ray behaviour. 
Note, however, that the discussion based on the X-ray lumi- 
nosity estimated from the simulations only gives the neces- 
sary condition for the quiescent state. 

In order to see whether our model satisfies a more strin- 
gent condition, we also have estimated the magnetospheric 
radius by using a general formula with the mass-accretion 
rate from our simulations and compared with the corotation 
radius for 4U 0115+63. The result was that our models with 
r < 1.2 are consistent with the observed X-ray behaviour. 
Other models with F = 5/3 and fc<0.6 were ruled out, be- 
cause in these models the system would show an X-ray out- 
burst after every periastron passage, which has never been 
observed. We have to admit, however, that the radius of our 
inner simulation boundary is about 60-100 times as large 
as the corotation radius of the neutron star in 4U 0115+63. 
To make the constraint on the system parameters more re- 
liable, we need numerical simulations with a much smaller 
inner boundary. 

It is important to note that our model is not only con- 
sistent with the quiescent state of 4U 0115+63. It can also 
give a natural explanation for a short series of Type I X- 
ray outbursts after occasional Type II X-ray outbursts. The 
mass-transfer rate from the Be disc in the X-ray outbursts is 
likely much higher than in the quiescent state. Suppose that 
the mass-transfer rate from the Be disc in 4U 0115+63 is 
temporarily increased by an order of magnitude. Then, the 
estimated X-ray luminosity in our model increases by an 
order of magnitude to a level of several x 10'^'' erg s~^ , a typ- 
ical X-ray luminosity of Type I X-ray outbursts. At the same 
time, the magnetosphric radius Rm decreases by a factor of 
about two, for which Rm is almost always smaller than the 
corotation radius, Rq, if fc ~ 0.5 (see Fig. llSf b'l'l. This means 
that the system is in the direct accretion regime. Thus, our 
model explains the Type I X-ray outbursts in 4U 0115+63 
as a temporal phenomena with an enhanced mass-transfer 
from the Be disc. 

Be/X-ray binaries are distributed over a wide range 
of orbital periods (16 d < Porb ^ 243 d) and eccentricities 
(e<0.9). It is interesting to consider here the effects of or- 
bital parameters. In this paper, we have seen that the main 
ingredients which determines the X-ray behaviour of a sys- 
tem are the viscous time-scale relative to the orbital period 
and the mass-transfer rate from the Be disc. Since the former 
is insensitive to the orbital period itself (see equation Q), 
the latter is most important if the viscosity parameter and 
the equation of state are similar in accretion discs of Be/X- 
ray binaries. 

In systems with ass ^0.1, the viscous time-scale in an 
outer part of the accretion disc is much longer than the 
orbital period unless the cooling is ineffective and the gas 
obeys an adiabatic equation of state. In these systems, per- 
sistent accretion discs are expected to form. Then, the mass- 
transfer rate from the Be disc determines whether the sys- 
tem shows Type I X-ray outbursts. For example, for a sys- 
tem with a moderate orbital eccentricity of e = 0.34 and 
a typical density of Be disc, the peak mass-transfer rate is 
~ 2 X 10~ ^'^ Mq yr"'^, for which the system stays in quies- 
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cence. It is obvious, however, that the situation should be 
different in systems with much higher rate of mass-transfer. 
A high mass-transfer rate is expected to result from a high 
orbital eccentricity and/or an intrinsically dense Be disc. 
Therefore, it is not strange that many systems which have 
regularly shown Type I X-ray outbursts have orbital eccen- 
tricities higher than that of 4U 0115-1-63 studied in this pa- 
per. In a subsequent paper, we will perform numerical sim- 
ulations of accretion in systems with high orbital eccentric- 
ities. 

Above we have assumed that the viscous time-scale is 
much longer than the orbital period in all Be/X-ray binaries. 
However, it is possible that an accretion disc has a relatively 
short viscous time-scale by some reason. One possibility is 
that the disc has a low cooling efficiency and becomes hot, 
particularly in an inner region. Our model 4 with F = 5/3 
gives an extreme example of such discs (see also Fig.jHJ. As 
shown in Fig. 1161 in a system with a short viscous time-scale, 
the mass-accretion rate has a high and sharp peak shortly 
after every periastron passage. The resemblance between 
the accretion rate profile in model 4 and the X-ray light 
curves of regular Ty pe I X-ray outburst s from systems such 
as EXO 2030-1-375 llWilson et alJ l2002l) implies that these 
are systems with relatively short viscous time-scales. 

Be/X-ray binaries are an ideal group of objects for 
studying the physics of accretion under the circumstances 
which have been studied little. Unlike close binaries, many 
Be/X-ray binaries have eccentric orbits. They can also be 
highly inclined systems. Such systems provide us a valuable 
oppotunity to study the effects of the periodically-changing 
tidal potential and mass-transfer rate and the inclination an- 
gle on the structure and evolution of the accretion flow. An- 
other important characteristic of Be/X-ray binaries is that 
they are systems with double circumstellar discs (the Be disc 
and the accretion disc). The interaction in such systems is 
threefold: the interactions between the neutron star and the 
Be disc, the Be star and the accretion disc and the Be disc 
and the accretion disc (mainly via the mass transfer) . Much 
more work is desirable both theoretically and observation- 
ally in order to understand this interesting and important 
group of objects. 
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